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Abstract 

We consider the problem of two capacitively coupled Josephson junction ar- 
rays made of ultrasmall junctions. Each one of the arrays can be in the semi- 
classical or quantum regimes, depending on their physical parameter values. 
The former case is dominated by a Cooper-pair superfluid while the quantum 
one is dominated by dynamic vortices leading to an insulating behavior. We 
first consider the limit when both arrays are in the semiclassical limit, and 
next the case when one array is quantum and the other semiclassical. We 
present WKB and Mean Field theory results for the critical temperature of 
each array when both are in the semiclassical limit. When one array is in the 
semiclassical regime and the other one in the quantum fluctuations dominated 
regimes, we derive a duality transformation between the charged and vortex 
dominated arrays that involve a gauge vector field, which is proportional to 
the site coupling capacitance between the arrays. The system considered here 
has been fabricated and we make some predictions as to possible experimen- 
tally measurable quantities that could be compared with theory. 

Key words: Superconductivity, macroscopic quantum phenomena, mesoscopic 
systems. 
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Personal note on Leo Kadanoff 



I am delighted to be able to add this contribution in honor of Leo Kadanoff's sixti- 
eth birthday. I have purposefully chosen a subject that involves duality transformations. 
Although in a completely different context, this goes back to several ideas I first learned 
from Leo many years ago, when I started my research career under his tutelage. Duality 
transformations continue to lead to new and interesting physics. Leo has been involved with 
duality transformations in several different contexts, making seminal contributions to their 
understanding and applications. 

I spent about three years with Leo, two as a his Ph.D. graduate student and one as his 
postdoctoral fellow. During this time I was able to see a real Master at work. Although 
we did not meet very often, every time we did (and he was always available to meet me), 
his insights and quickness left a deep impression on me. Everybody has his/her own style 
of doing research, and while I could not claim to do research as Leo does, his approach, 
style and his depth made a lasting impression on my way of looking at physics and physics 
problems addressing their detailed quantitative solutions. 

I. INTRODUCTION 

The subject treated in this paper relates to interesting quantum properties of Josephson 
junction arrays (JJA) made with ultrasmall junctions. Layered Josephson junction arrays 
have been the source of many theoretical and experimental studies in the last few years 
|IJ. Recent advances in submicrometer technology have made it possible to fabricate rela- 
tively large arrays of ultrasmall superconductor-insulating-superconductor (SIS) Josephson 
junctions The areas of these junctions can vary from a few microns to submicron 

sizes. Under these circumstances the long range phase coherent properties of the JJA de- 
pend crucially on the interplay between the Josephson energy, Ej, and charging energy, Eq- 
Detailed experiments have been carried out, for example at Delft, that have produced a 
phase diagram of temperature vs the quantum parameter a = ^ @. We have calculated 
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the a vs temperature phase diagram and we made |||7[ a direct successful comparisons to 
the experimental results . Our results were obtained using a WKB-renormalization-group 
approach, plus a variational and quantum Monte Carlo (QMC) calculations ||. An im- 
portant QMC result is that there appears to be a low temperature QUantum fluctuation 
Induced Transition (QUIT) in this system 

For the most part the experimental systems have been two-dimensional, but prototype 
quasi-three-dimensional samples have also been fabricated ||. In this paper quasi- means two 
layers of JJA capacitively coupled at each lattice site. There are two dominant contributions 
to the charging energy in the type of junctions fabricated ||; one due to the addition of 
a Cooper pair charge to a superconducting island given by E Cs = jt- with C s the self 
capacitance and e the electronic charge, and the charging energy necessary to transfer a 
Cooper charge from one island to its nearest neighbor, given by Ec m = |f-, with C m the 
mutual capacitance. In the Delft experiments, C s ~ 3 x 10~ 17 F, and C m ~ 1 x 10 _15 F, 
which mucins tlicit C-m can be two orders of magnitude larger than C s . 

The phase diagram for one layered JJA has the following general characteristics. At 
low temperatures, for small a there is a superconducting phase in which the Cooper pair 
charges are delocalized, while for large a the system has delocalized vortices and the array 
is an insulator. There is a phase boundary that separates the superconducting to insulating 
regions. Now assume that we have two layered JJA capacitively coupled at each lattice site. 
This configuration is potentially quite interesting since, as mentioned above, each array 
can be in one of two extreme limits; one Cooper charge dominated and the other vortex 
dominated. Each array is now described by its quantum ratio «j = ^r 1 (i = 1,2). We can 
then imagine to have the two arrays in four possible configurations. When «i < 1, the i-th 
array is dominated by localized vortex excitations, Vi, while the Copper pair excess charge 
excitations, Qi, are in a superfluid or superconducting state. In the «j ^> 1 regime the array 
has the Q,'s localized in an insulating state while the V^'s are delocalized. We can also have 
the extreme cases when both arrays are semiclassical or in the quantum regimes. 

In this paper we shall consider first the case where both arrays are in the semiclassical 
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regime. We obtain an effective partition function that allows to calculate the change of 
each array critical temperatures as a function. This type of analysis was quite informative 
in the one-array problem. Next we move to consider the general case where we derive a 
Hamiltonian that is valid for all parameter regimes. This Hamiltonian is quite complex and 
difficult to analyze in full. Instead we consider the interesting case when one of the arrays 
is charge and the other vortex dominated. After a series of transformation we arrive at an 
effective Hamiltonian that exhibits interesting duality properties. We analyze the imaginary 
time dynamics of this Hamiltonian in the case where only one vortex in one array and one 
charge in the other are considered. Here we show that these tow excitations interact via a 
gauge-like interaction proportional to the interaction capacitance between the arrays. 

The outline of the paper is the following, in Section III the case where both arrays are 
in the semiclassical limit since in that case we can analyze in some detail the changes in 
the critical temperature of both arrays independently via a variational Mean Field analysis. 
In Section IV we consider the interesting limit when one array is semiclassical and the 
other quantum. In this case we can derive an effective action for the problem that allows 
a general analysis of the interaction of a vortex in one array with a charge in the other 
via an interaction term that has the form of a minimal gauge coupling proportional to the 
interaction capacitance. The other extreme case where both arrays are quantum is not 
considered here since that case is harder to analyze. We conclude the paper with some 
conclusions in Section V. 

In this paper I present some results of work done in collaboration with C. Rojas (related 
details can be found in || and in his Ph.D. thesis 0). 

II. THE MODEL 

In this section we define the model that describes a quasi-three-dimensional array com- 
posed of two JJA layers coupled at each site by an ultrasmall capacitor. In our analysis 
we will have in mind the prototypical samples fabricated at Delft M. In these samples the 
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size for each layer were L x = 230, and L y = 60. The typical parameter for the intra-layer 
mutual capacitance was C m ~ 2.3JF, and the interplane local interaction capacitance was 
Cint ~ 0.6/F. The model consists of two planar arrays stacked on top of each other. The 
intra-array interaction between the superconducting islands in each array contains an elec- 
trostatic and a Josephson coupling. The two planes are only capacitively coupled. The 
model Hamiltonian is given by 

H = \ E E V M (r>„(r) + F^}) + F 2 ({0 2 }), (1) 

where the index /i = 1,2 labels the two arrays. The operators M and n M satisfy the 
commutation relations [n^ri), ^'{r^)] = —iS^^d^, /J,'. The functions F M are the Josephson 
interaction terms, 

F,m) = Ef E (l-cos(Mn)-Mm. (2) 

<n,r2> 

Here is the Josephson energy coupling constant, for the junctions in array fi. V M (r*) 
is the electrostatic potential felt by the charges contained in the superconducting island 
located at r in array fi. This potential is produced by all the other charges in both arrays, 
and it is obtained from the discrete Poisson equation 

q m(f) = Cg) E [Vi(r) - Vi(f + «)] + Ci^ViCf) + C int [V 2 (f) - Vi(f)] • (3) 

u 

Here the w summation is over nearest neighbors. In a uniform square lattice u = {±x, ±y}. 
The complementary equation for n 2 (r) is obtained by interchanging (1 <-> 2). In Eq.(^), 
C s ^) is the self capacitance of a superconducting island in array \i. The previous equation 
can be written in the compact form, 

Q n^ifx) = E E C^(ri, f 2 )V„(f 2 ), (4) 
where the capacitance supermatrix C^ u is made of four blocks labeled by /i, v = 1, 2, 
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(5) 



(CM + *C<J*> + C int ) , if /i = z/ and n = f 2 , 

—C$ , if /i = z/ and r*i = r 2 ± u, 

-Qnt, if n^v and Fi = r 2 , 

0, otherwise. 

The diagonal blocks of this matrix are the intra-array capacitance matrices. We will use the 
following notation for them 



C M — C^. 



(6) 



The off-diagonal parts of the supermatrix are given by block matrices proportional to 
the identity matrix — C'j I1 tIzv,jv'j with N the linear size of the arrays. The inverse matrix of 
Eq.(||), C, is obtained from solving the equation 



EE C /^(^i>^2)C yjP (f2,f 3 ) = 5^ p 6? u 
The explicit components of the inverse matrix are given by 



r-.i ■ 



(7) 



Ci,i - C x 1 I - cf nt c 2 1 C 1 1 



'2,2 



cr 1 



i-c^c^c, 1 



Ci,2 — C- m tC 1 1 C 2 I — C( nt C 1 i C 2 



C 2) i — Ci nt C 2 1 C 1 I — C\ at Q 2 i C 1 



(8) 
(9) 
(10) 

(11) 



For a uniform array an explicit expression for this matrix can be found using a Fourier 
representation of the matrices. 

Finally, using Eqs. (jl]), @, and (|7|), we can write down the model Hamiltonian studied 
in this paper as 

H = y E E ^(n)C^(n, F 2 )n,(r 2 ) + F 1 ({} 1 }) + F 2 ({0 2 }). (12) 
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III. SEMICLASSICAL CAPACITIVELY COUPLED ARRAYS 



In this section we begin our analysis of the quasi-three dimensional JJA model defined 
by Eq.flT2|), when the parameters in both arrays are in the small a semiclassical regime. 
Here we estimate the change in the critical temperature, T c , for each array as a function of 
Cint- We use a WKB semiclassical expansion valid for small a values, as we did in the single 
array problem ||||. In this limit we first obtain an effective action that we then analyze 
within a variational Mean Field theory (MFT) approach. In the prototypical Delft samples, 
they had small ratio values for (Ci nt /C s ). To have a consistent semiclassical expansion in 
Cint, as follows from looking at Eqs.(||) -([□]), we nee d to carry out the expansion at least 
to second order in C, which is equivalent to doing the expansion up to second order in q 2 . 
After a long, but direct, calculation we obtain the effective semiclassical partition function, 



j sc 



exp{-S cfi [(f)]}, (13) 



with the effective semiclassical action up to second order given by 

= /?£ ji - ^ [c,A\o\) - c^(\u\) 
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C^(|0|)-C^(|«|) 



+ 720 - d<pi C - wj- (14) 

As a check, we note that we recover the one-layer results if we keep terms up to first order 
in q 2 . We find that the superconducting state is stable at low temperatures, because the 
second-order contribution in q 2 has a negative sign compared to the first-order contribution 
to the effective action. Notice that higher orders in q 2 coincide with higher orders in f3. This 
conclusion does not, by itself, eliminate the possibility of having a QUIT in the two layered 
problem. 

Eq.(|i"4j) does not have a simple form as in the one JJA layer first-order expansion in 
q 2 . The third term in Eq.(|T4|) has nonlocal interactions, which makes a direct calculation 
of the critical temperature in general complicated. To estimate the change in the critical 



temperatures as a function of C- m t, we have performed a couple of distinct MFT variational 
calculations for the partition function given in Eq. ([H) . 

For the general variational calculation first we split the effective action into two parts 

s^faM = s^[<Pi] + 5j 2) [0 2 ] + - sf - 4 2) N) , (15) 

which gives the exact semiclassical expansion for the partition function 

Z = Z« 4 2) (exp {-5W[^l, 02] + + 4 2) [0 2 ]}) o , (16) 

Z^ = Jd^ exp{-5 [^]}. (17) 

The average <> is defined as 

( A )o = ^n^® ffyidfaAtyufa) ex P {-Sf (18) 
A) z o 

We now use the variational inequality {exp{A}) < exp{(A)} to write 

Z<Z^ 4 2) exp{(-S eS [ ( j )l A2]+S^ 1 }+S^[ ( p 2 }) o }. (19) 

The corresponding variational free energy is then 

PF var = -InZ^ - InZf - (-S eS [<f> 1} 2 ] + 5« + S^[<h]) Q ■ (20) 

We next need to specify the functions Sq . They can be any general functions, but here 
we restrict them to be functions of only one phase variable for each layer. These functions 
must be chosen so that we can carry out some or all of the integrations in Eq. (|2~0"D . We 
introduce variational parameters in the trial actions that are determined by requiring that 
they minimize (3F wax . We have used two different choices for S^ 1 . 

A. First variational calculation 

The first choice decouples all the phases in both arrays 

# ) foj=-7,£cos(^(f0), (21) 

r 
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with the two 7 M 's the variational parameters [|10| . The advantage of this form for the trial 



action is that all the integrals can be analytically computed or they can be expressed as 
simpler one-dimensional integrals. For the classical 2-D XY model it is known that this 



variational choice grossly overestimate the critical temperature |]I0 |. This approximation, 
nonetheless, gives good qualitative results for the critical temperatures for both JJA layers. 
Here we are mostly interested in general trends so, for simplicity, we only study the case when 
both JJA layers have the same parameter values. After a lengthy but direct calculation, the 
critical temperature equation for one of arrays is obtained from solving the equation 

D A x 4 + D 3 x 3 + D 2 x 2 + D 1 x + D = 0, (22) 

were we found 

x = (3 c Ej, a = e 2 /(2C s Ej), D = -1, D, = 2, (23) 



D 2 = -(4/3) atfx, D 3 = a 2 K 2 /9, D 4 = W{K 2 /2 - K 3 )/45. (24) 

Ki,K 2 , and K 3 are complicated functions of the ratios (C- mt /C b ) and (C/C s ). They can 
be computed in general in terms of summations over Fourier modes. Here we are only 
interested in a general trend of T c as a function of C- mt , and we further simplify the problem 
by considering the self-capacitive limit, i.e. C m = 0. In this case, the K functions can be 
fully computed giving 

K ^ = ( \lfc!!c) > K * = 28K l and K 3 = 8K 2 . (25) 

We can now interpret the effect of the inter-plane capacitance as a rescaling of the single 
layer quantum parameter a, namely 

a " = (i + kwcJ a - (26) 

The first conclusion we draw from this result is that, in the semiclassical approximation, the 
inter-plane capacitance makes the system less quantum mechanical. The critical temperature 
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increases as (Ci nt /C s ) increases up to an asymptotic plateau. To further check this result we 
have also performed quantum Monte Carlo calculations that confirmed our analytic results. 

Note that here we only presented the case where both arrays are equal and the mutual 
capacitance is zero. We have done so due to the simplicity of the analytic results. It is not 
much harder, within this approximation, to numerically calculate T c for the more general 
cases. The general result leads to the same qualitative conclusion; The increase in the 
inter-plane capacitance raises the critical temperature for both arrays. 

From Eq.([26l) we can find T c up to second order in the effective quantum parameter a e s 
giving 

(k B T c \ (k B T^\ /2\ /189\ , 3 , 

{-et) = [it) - (3) a " + (wo) + 0(Q '* ) ' (27) 

where the variational result for the classical 2-D XY model is (ksT^ ' /Ej) = 2. It is evident 
that this approximation overestimated the critical temperature. Surprisingly, it gives a very 
good estimate of the first two correction values for the one array problem M. 



B. Second variational calculation 

In the previous subsection we used Eq.([H]) to perform the MFT variational calculation. 
From the form of Eq.(^) it is clear that it is not necessary to use a trial action that decouples 
all the phases in both arrays. We now want to use a better trial action that gives better 2-D 
XY model classical results. This has the advantage that the classical limit for each array is 
by construction exact. The disadvantage of this choice is that we need to evaluate a nonlocal 
average over the 2-D XY model classical Hamiltonian, for which we have to approximately 
evaluate an infinite lattice sum. Again, we restrict the calculation to the C m = limit, 
although we could do the numerical calculation for the full model. 

The starting point of this scheme is to use the trial action 

3^^] = -^ E cos(^(ri)-^(r 2 )), (28) 

<n,r2> 
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with fti, and fa the variational parameters needed to calculate the variational free energy in 
Eq. (pop from Eq. (|i4|) . We need to evaluate the following averages 

(F fl ) XY = E { f £ ^-(cos(Uri)-Mr2))) XY ), (29) 



<ri,r 2 > 

^f) ~ (cos (^(rO - M (f 2 )) cos (0„(f 3 ) - <^(r 4 ))) xy 

(30) 



The average in Eq. (^) is simple, since we only need to perform the average when T\ and fj 
are nearest neighbors. Here we are interested in periodic and symmetric arrays. It is then 
convenient to define and evaluate the following short range correlation function 

go(j3) = (cos(0(f) - 0(f + u)) xy , (31) 

with (3 is the inverse temperature of the classical 2-D XY model. This function can not be 
calculated exactly, but it can be evaluated using a Monte Carlo calculation or by matching 
a low to a high-temperature expansion for the classical 2-D XY model. 

The average in Eq.(pO|) is more complicated. First, if \i ^ v, due to the intra-plane 
independence, the average can be reduced to finding two values of go(fl). If \i = v, the 
problem is intractable for general f\ and r 3 . This is exactly the problem we would encounter 
if we want to make a calculation for a general full capacitance matrix. On the other hand 
if we choose C m = 0, the problem is simplified. In this case, only go(/3) and the following 
functions 

9l {(3) = ([cos (0(f) -0(f+x))] 2 ) xy , (32) 

9ziP) = ( cos (0(f) - 0(f + x)) cos (0(f) - 0(f + y)) ) xy , (33) 

o [ 2{dg x /d(3) + (dg 2 /dj3) \ 
9 ^ ) = 2 { dg^/dp /' (34) 

need to be known. These function can be calculated using a Monte Carlo calculation for 

the classical 2-D XY model. Using all these functions, we can evaluate the variational free 

energy F w (/3i, ^2). The parameters Pi and are found by imposing the condition that 

they minimize the free energy. This condition reads 
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8F BF 

^T = 0, ^ = 0- (35) 
df3i d(3 2 

The explicit equation for fa obtained from Eq.(|35"D is 

a - (bf^U , iP E ?)\ nA , 4(M 1} ) 4 f (1) 2 

Pi - IP^J JS 1 3 "off "I «eff ) H ^ 1 "eflf ^3 (Pi J + 

2 



If I £ n (!) (2) ( Cint/CP \ ( C int / C s (2) t i 



A similar equation can be written down for [3 2 . In writing this equation we have used the 
following definition for the effective quantum parameter for the arrays 



ri) _ e 2 ( 1 + C int /CP . 

"" !! ■~2c*M ) Vi + c int /a (1) + a nt /d 2) ^ 1 j 



We are therefore left with the following set of self-consistent equations for the two variational 
parameters 

Pi = Gi(/3,/3i,/3 2 ), P2 = G 2 (P,Pi,P2). (38) 

The functions Gi and G 2 can be identified from Eq. (|36|). We set /3 M = /3 ( f y as the condition 
to find the critical temperature of the array //, i.e. we identify the value of /3 M with an 
effective inverse temperature for the array [i. To find the critical temperature for array 1 we 
solve the following set of equations 

Pc Y = GM 1 \pr,P2), h =G 2 {$\ff Y ,fr). (39) 
These equations are not difficult to solve. First, given a value of /3 we first solve for /3 2 using 



Eq . (p9|) . Note that this equation can be written in the following form 

fc = A + B g 3 (p 2 ). (40) 

We can show from Eq.([3(J) that the quantities A and B in this equation are not functions 
of Moreover, solving this equation using a fixed point search method gives {3 2 (f3^), that 
can then be introduced into Eq. ([39|). We are now left with a one variable equation. 
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The case when both arrays have the same parameters is again easier to solve, since the 
transcendental Eqs. (|39| ) reduce to just one polynomial equation 

RiX 4 + R 3 x 3 + R 2 x 2 + R x x + #0 = 0, (41) 

where /3« = and 

2 1 

Ro = —f3 c , #1 = 1, #2 = -gtteff, R3 = Yg^eff' ^ 42 ) 



p 4 2 



3.6 + 9.6 



Cint/ Cs 



(43) 



, 1 + Cint/ Cs / 

We have used Eq.(p6|) to define a e s and found the value of gz{fic Y ) using Monte Carlo 
calculations to be 

tf Y ^ 1.1186, ^ 3 (/5 c Xy ) « 3.6 . (44) 

The result of this calculation gives the same result as in Eq.(^) for the first order correction 
in a, but now the a = limit gives the correct classical value for T c /kBEj = 1/ (3^ Y . Here 
the general result is qualitatively the same; an increase in the interaction capacitance C- m t 
results in a decrease of the effective quantum parameter. 



IV. DUALITY IN TWO CAPACITIVELY COUPLED JJA 

In the previous section we studied the specific changes in the individual critical temper- 
ature of two JJA, when both were in the semiclassical parameter regime. In this section we 
analytically consider the interesting case where one array is in the semiclassical regime and 
the other in the quantum one. In the semiclassical array, vortices are localized while the 
Cooper pair charge fluctuations are mobile in the superfluid phase. This happens for small 
a values. In the quantum JJA regime, large a, the vortices are mobile while the charges 
are localized. It is not possible to have both vortices and charges simultaneously localized 
or mobile since this is forbidden by the Heisenberg uncertainty principle (which has been 
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shown at work in recent experiments in a small Josephson array system [0). The inter- 
action between vortices and charges has a minimal coupling form, with constant strength 
and it is sharply localized, i.e. a vortex and a Cooper pair only interact if they are located 
at the same point in the array. By considering two arrays with one vortex dominated and 
the other charge dominated we can have them both interacting via the coupling capacitance 
between the arrays. A related analysis of the coupled array system was considered in Ref. 
15fl by us and an alternative and complementary analysis was also presented in Ref. . 



In this section we will carry out the two array analysis extending techniques developed 
for the study of one array ||12|| . Here we shall consider first the one-array component of the 
Hamiltonian given in Eq.([E|) and later its Villain approximation | 13| . We briefly mention 



the one JJA calculational approach, since the extension to the two-array case follows from 
this analysis. This is true only because the arrays are electrostatically coupled. 
The Hamiltonian for one array reads 

H=^J^h{f 1 )C-\r l) r 2 )h{r 2 )+E J £ [l - cos (fe) - 0(f 2 )) ], (45) 

ri,r*2 <fi,T2> 

where C(fi,f2) is the one-array capacitance. The partition function for this Hamiltonian 
can be written in the path integral form 



nn e / -^-p -y d - 



{n(r,r)} 



q 2 

-^]n(r,ri)C 1 (r 1 ,r 2 )n(T,r 2 ) + 

fx ,?2 



(46) 



Here we denote 0ri,r 2 ( r ) = < f ) { T > n.) ~ ( t ) { T -, ^2)1 an d write the imaginary time summation as an 
integral. To integrate over the 0's, we need to introduce an additional set of variables. This 
is done by writing the Boltzmann factor as a Fourier series using the Poisson summation 



representation JL3 



exp{-A(l -cos^)} = /m(A) exp{z^m}, (47) 

m=— 00 

/•2tt dip 

/m(A) = / — exp{-A(l -cos-0)} exp{-i-0m} = exp{-A} I m (\). (48) 

JO 2,71 
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Here I m (X) is the modified Bessel function. The integral in Eq.fl48|) gives a convenient way 
to extract the asymptotics for the small and large A. For small A, the Taylor series expansion 
of 7 m (A) around A = gives us the leading expansion terms of f m (X). For large A, a steepest 
descent calculation yields the leading term. In these two asymptotic limits we have 



/m(A) 



(A/2) m /m!, 



if A < 1, 



(49) 



exp{-m 2 /2A}/V2^A, if A > 1. 
If we use the A ^> 1 result in Eq.(f|7|), we can write 

1 



exp {—A (1 — cos ■?/>)} 



exp{— m 2 /2A + iiprn]. 



(50) 



Using Eq. (|47D in Eq.(^) we have new summation link variables m v {r, r) between the nodes 
of the lattice. This is repeated for each of the imaginary time planes. The subindex v is a 
vector that denotes the orientation of the links so that we can write 



exp 



{ - eEj [l - cos (0(r, r + v)- (j)(r, r)) } 



/ 1 



\2iieE, 



x 



X 



6XP { 2eEj 



m v (r, rf + im u (r, r)A u <p(T, r) 



(51) 



Here we have discretized the imaginary time interval with e = /3h/L T , and the lattice 
derivative is A u f(r) = f(r+ v) — j{f). After the integrations over the 0's the result can be 
written in terms of the integer n and m variables as 



Z ~ E Eexp -E 

{n(r,r)} {m{r,r)} I t 



(52) 



After the integration over the 0's we get a set of constraints over the n and m values. These 
constraints can be written as discrete continuity equations that are satisfied at each node of 
the array, 



A ■ m(r, r) + A T n(r, r) = 0. 



(53) 
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These constraint equations can be solved in several ways. For example, the pair (n, to) can 
be expressed in terms of a three- vector /C = (n, m) ]I1J], so that Eg. (|53|) becomes 

A„/C" = 0, (54) 

with the discrete gradient, A u = (A r , A). From Eq.(Q) K can be expressed as the curl of a 
gauge field, i.e. /C M = e pvp /S. v A pi with e pvp the usual fully antisymmetric tensor. Substituting 
this result into Eq.(^) we get an effective action over the gauge field A, which resolves the 
constraints over the summations. 



Our solution to the constraint in Eq . (p3|) is different from the one used in Ref. |T2 
where they wanted to preserve the n variables. The solution to our constraint equation will 
have a particular solution plus an homogeneous solution. Note that Eq.(^3|), written in this 
form, resembles one of Maxwell's equation that connect the divergence of the electric field 
to the charge density. The particular solution to this equation contains the gradient of a 
line integral, that can then be solved using a discrete line integral operator. The solution 
can be formally written as 

m^r, r) = -e"(e • A) -1 A T n(r, f) + e^A u A(r, f), (55) 

with A(t, f) another integer gauge field. To obtain the partition function we need to perform 
summations over this field. The first term in the last equation represents a discrete line 
integral, which is better calculated in a Fourier representation. After substituting Eq.(|55"D 
into Eq.(|52D we obtain an expression for an effective action as functions of the n's and the 
A's. The partition function is now obtained by summing over these variables, but now 
without any constraint. We perform the summation over the A's again with the help of 
the Poisson summation formula. After introducing a new set of integer v variables, the 
summation over the A's can now be done since the integrals left to calculate are Gaussian 
and the A's are unconstrained fields. The final result is 

^ = EE ex P - S cS (n,v) , (56) 
M M L J 
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S eS (n,v) = E 



q 2 e 
~2i7 



n(r, ri)C(fi, f 2 )n(r, r 2 ) + 7ieEjv(r, fi)G(n, r 2 )v(r, r 2 ) + 



+i n(r, fi)0(fi,r 2 )A T u(r, f 2 ) + 



AixeEj 



A T n(r, Fi)G(Fi, f 2 )A T n(r, r 2 ) 



(57) 



where we defined 



G(ri, r 2 ) ~ In Ifi — r 2 |, anc/ ©(r*i, r 2 ) ~ arctan 



(58) 



.Xi - x 2 , 

Eq. ( p7[ ) is an effective action for two coupled imaginary time Coulomb gases. This equation 
is valid for all parameter ranges. When Ej is large, the last term in this equation is small and 
the time derivatives of the charges are soft with the ra's having strong fluctuations. In this 
limit the v's dominate and the n's are not well defined. We will call this a vortex-dominated 
regime. When Ej is small, the last term in the effective action is large and it makes the time 
derivatives of the charges well defined. In this limit the v's are not well defined and the state 
is charge-dominated, with the charges described by an effective continuous Gaussian model. 
After integrating the continuous variables we obtain an effective action for the vortex integer 
conjugate variables. 

One important aspect of the one-array derivation of Eq. (p7[) is that it did not involve 
the charging energy part. This means that when we carry out the two-array calculations 
we only need to see that its effective action can be written down immediately from just 
repeating the one array calculation; we only need to add the essential extra charging energy 
term that couples the two arrays. The two-array equivalent equation to Eq.([57]) is then 



s ef Kn«yv {2) y 2) ) = E 



— n (1) (r, n) Ci i(fi, f 2 )n (1) (r, r 2 ) + 

271 

+7re£i 1 V 1) (r, ri)G(f 1; f 2 )v {1 \r, f 2 ) + 

+i n (1) (r,fi)e(fi,f 2 )A T w (1) (r,f 2 ) + 
1 



+ 



+ E 



47re£ 
2vr 



(i) 



A T n {1 \r, ri)G(fi, f 2 ) A T n {1) (r, f 2 ) 



+ 



n (2j (r, rjjCi^ri, r 2 )n (2) (r, r 2 ) + 
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+neEf ] v^(T, n)G(fi, f 2 y 2) (r, r 2 ) + 



(2), 



+i n ( - 2 \r,f 1 )&(r 1 ,f 2 )A T v {2) (r,r 2 ) + 
1 



+ 



A-neE 



A r n( 2 )(r,f 1 )G(f 1 ,f 2 )Ay 2 )(r, 



r 2 



+ 



+ E 



— n (1) (r, fi)Ci j2 (ri, r 2 )n (2j (r, r 2 



(2), 



7T 



(59) 



where we used the definition of C given in Eq. (|7|). 

In this section we shall consider the interesting case where one of the arrays is in the 
semiclassical regime (the vortex-dominated state) and the other is in the full quantum regime 
(i.e. the charge-dominated state). We want in particular to study the interaction between 
vortices in one array and charges in the other. We assume, as in experiment, that the 
arrays are dominated by the mutual capacitance between nearest neighbors. We take array 
1 vortex-dominated and array 2 charge-dominated, i.e. 



(60) 



We start by performing a vortex integration in array 2. This can be done using the Poisson 
summation formula to write 



{v( 2 )} 



s/n d *( T .*o«p[-5«($) 

2m^2^(r, r)P(r, r) 



{P} r,r 

x exp 



(61) 



(2) (2) 

We can neglect the P^O terms when Eq 3> Ej , since they are exponentially small. Note 
that this integration only involves vortices in array 2, and it can be done without affecting 
the variables in array 1. When we only consider the P = terms, there is a change in the 



Josephson coupling constant given by Ef E K j'/2 |g . In this parameter limit wc sec 
that the only modification is in the new coupling constant in Eq. (|5TS|) . 

We can also use the Poisson summation formula for the integration of the charges in 
array 1. After the integration, the action for the charges can be written as 



,(2) 



T,T f s ri ,7*2 



n 



(1). 



r, ri)M(r, r'; f h f 2 )n (1) (r', f 2 ) + £ J(r, r>«(r, 



(62) 



where the M operator and effective current J are 



a 2 e ~ 1 

27reE {1) 



r 2 L 



q 2 e~ 



ie(f 1; r 2 )A rU W(r,f 2 ) + ^C 1>2 (f 1; r 2 )n( 2 V,r 2 ) 

7T 



(63) 
(64) 



After integrating the charges in array 1 and the vortices in array 2 we are left with the 
following expression for the effective partition function 



Z = £ £ exp[-S cS (v^\n^) 

{vW} {n( 2 )} 

where the effective action for vortices in array 1 and charges in array 2 is given by 



S cS [v {1 \n^]= Y, 



T,r\,T2 



' A r n( 2 )(r,f 1 )G(f 1 ,f 2 )Ay 2 )(r, 



r 2j 



+ 



+ E 



T,T,ri,T2 L 



-n (2) (r, fi)G„(r, r'; f 1; f 2 )n (2) (T', f 2 ) + 



+i n (2) (r,fi)e(r, r'; r b r 2 )A r w (1) (r', r 2 ) + 



f-^A^^ 1 ) (r, fi)G v (r, r'- f u f 2 )A T v^{r', 
2q z e 



r 2 ) 



The effective interaction potentials are defined by 



G n (r,r';fi,f 2 ) = — 

7T 



C 2 , 2 (ri,r 2 )5 TiT / - 

E Ci, 2 (ri, ?3)M-\t, t'; f 3 , f 4 )C 1)2 (r 4 , r 2 ) 

e(-.-':n.r2) E e(r : 2 ,^)M- 1 (r,r / ;^,r 4 )C(r : 4,r : 2 ), 

G v (r,r';ri,r 2 ) = £ 0(ri, r 3 )M- 1 (r, r'; f 3 , f 4 )0(r 4 , f 2 ). 



g 2 e 

7T 

g 2 e 

7T 



(65) 



(66) 



(67) 

(68) 
(69) 



Notice that the time nonlocality of these kernels comes from the second term in Eq.(|63|). 
To gain some physical understanding of these complicated equations we will next discuss a 
simplification when the nonlocal term in Eq.(|63"D is small. 
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A. Vortex-charge capacitive gauge-like coupling 



In Eq.(^6|) we have an effective interaction between vortices in array 1 and charges in 
array 2. We consider the dynamics of just one charge and one vortex in each array, in a 
similar way as was done for one array in Ref. [O]. Lets assume that the vortex and the 
charge move along the imaginary time- dependent trajectories R{t), and X{r) respectively. 
In this case the vortex and charge space-time distributions can be described by 

v^(r,f) = 5[f-R{T)]; n {1) {r,r) = S[f-X(r)}. (70) 

After taking the time derivative of v^ 1 ' we find 

A t vW(t,?) = A T S[n - R(r)] = -£A M <f[ri - R(t)] A T R(r). (71) 

The right hand side in this equation relates the time derivative to a summation over space 
derivatives. We can next rewrite the interaction term in the effective action in the following 

way 

rP % rP % ~ / - \ dR (r f ) 

S mt = -i j dr J o dr' £ A,0 (r, r'; X(r), R(r')) (72) 

This expression has a similar form to the typical minimal gauge coupling in electrodynamics. 
Following this analogy we can define the corresponding vector potential 

A(t, r) = ^ drAQ (r, r'; r, X(r)) . (73) 

Using this definition Eq.(^) can be rewritten as 

S int = -i J^dr' A(t',R{t')) ■ (74) 

Here we have chosen to view the vortex-charge interaction in the representation where the 
vortex moves under the influence of the charge gauge-like field A. This view is equivalent to 
the representation where the charge moves in the gauge-like field produced by the vortex. 
We have a vector field, and we can find its corresponding effective "magnetic field" 
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B(r,r) = A x A(r,r), 

= -- j dT 'J2 (AxA@(X(r'),r 3 ))M-\T,T , -,r 3 ,r 4 )C 1<2 (r i ,r). 



(75) 



The solution kernel for a point vortex at the origin satisfies the equation 

A x A0(ri,f 2 ) = 2vr 5? lf2 k, (76) 



b(r,f) = - q — / drXM-^r.rUfrVDCi^ri.r) k. (77) 
it Jo 



from which we can get our final expression for the effective magnetic field 

dr' 

r'l 

Up to now we have that the effective action in Eq. (|66|) , the effective gauge vector potential, 
and its magnetic field interaction are nonlocal in time, due to the nonlocality in time of the 
second term in Eq. fl6"3|). In the limit {f3E { 'P){f3E^ ; '^) 3> 1, the second term in Eq.(|63|) is 
negligible, and therefore we can write 

M-V.iVi.f?,) « ^-C^Kfx,^) 5 T y. (78) 

Using Eq.(D) we can write this equation in terms of the intra-array capacitance matrix 

C 1 ;i = [C 1 -C? t C 2 - 1 ]. (79) 

We note that the array is periodic and symmetric, so that all the commuting matrix operators 
can be diagonalized using plane waves. These facts produce important simplifications in the 
rest of the interaction kernels, giving the results 

e « -C int 0CJ 1 5 T y, (80) 



G n w q - Co" 1 5ry, (81) 



2 



7T 



G v « 0[C X - C int C 2 : ]0 5 T y. (82) 

Eq.(^Tj) is particularly significant, since it implies that, within this approximation, the in- 
teractions among charges in array 2 do not depend on the presence of array 1. This is 
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a counter-intuitive result, because we would expect that a virtual photon excited from an 
island in array 2 and absorbed in another island in the same array would have contributions 
from bounce interactions with array 1. What happens is that after adding all the contri- 
butions from these bounces, the net result (within the approximation leading to Eq.([7^)) is 
a cancellation of the contributions arising from array 1. Finally, we can write the effective 
magnetic field as 

b(^^-C int C^(f,X(T))k (83) 

This result implies that, if we have a charge at X(t) and a vortex at R(t), the vortex will 
feel an effective magnetic field produced by the charge of magnitude — C^C^ (X(t), R(r)). 
This situation is reminiscent of the vortex-charge bound states extensively studied in the 
fractional quantum Hall effect problem, and it may very well be that this system may serve 
as an experimental prototype for those types of problems. 

Our discussion here has concentrated on deriving and analyzing convenient partition 
function expressions that one can also use in quantum Monte Carlo simulations. We have 
done some work in this direction, but we must say that the problem is still highly non trivial 
because of the form of the kernels in the effective action. However, we expect to further 
unravel interesting physics for this problem in future. 

V. CONCLUSIONS 

In this paper we have introduced and presented results for a model of two capaci- 
tively coupled quantum Josephson junction arrays. This is a difficult problem but one that 
promises to lead to interesting new physics. We have first derived a semiclassical expression 
for an effective Hamiltonian, that allowed us to study the change in the critical temperature 
for each array. We used two types of variational actions that permitted the evaluation of 
critical temperature shifts as a function of the inter-layer coupling capacitance. The main 
qualitative result is that an increase in the interaction capacitance increases phase coherence 
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in the arrays. Next we considered the interesting case where one array is quantum phase 
dominated and the other Cooper pair charge dominated. Here we extended the one-array 



work of Fazio et al. |12| to the capacitively coupled two-array problem. We wrote an effec- 
tive action in terms of four interacting imaginary time Coulomb-like gases, and derived an 
effective Hamiltonian for the coupled system. The effective Hamiltonian is dually symmetric 
between charge and vorticity in form but with complicated kernels. In the simplified case 
where one array has one vortex and the other one charge, we showed that their interaction 
has a minimal gauge-like coupling. This interaction is, however, nonlocal in the gauge field. 

Finally, this type of system holds the promise to lead to a variety of novel experimentally 
observable macroscopic quantum phenomena. In particular, the vortex-charge interaction 
discussed at the end of this paper deals with the interplay of quantum-classical effects, and 
may lead to possible fractional statistics analogies to the fractional quantum Hall effect. 
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